Let us set some global options for all code chunks in this
document.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = TRUE,
# Disable warnings printed by R code chunks
warning = TRUE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)
# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
## Loading required package: Matrix
## This is INLA_25.04.16 built 2025-04-16 08:05:23 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - List available models/likelihoods/etc with inla.list.models()
## - Use inla.doc(<NAME>) to access documentation
## - Consider upgrading R-INLA to testing[25.04.29] or stable[24.12.11] (require R-4.5)
## Loading required package: fmesher
## This is rSPDE 2.5.1
## - See https://davidbolin.github.io/rSPDE for vignettes and manuals.
## This is MetricGraph 1.4.1
## - See https://davidbolin.github.io/MetricGraph for vignettes and manuals.
##
## Attaching package: 'MetricGraph'
## The following object is masked from 'package:stats':
##
## filter
library(grateful)
library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
We want to solve the fractional diffusion equation \[\begin{equation}
\label{eq:maineq}
\partial_t u+(\kappa^2-\Delta_\Gamma)^{\frac{\alpha}{2}} u=f \text {
on } \Gamma \times(0, T), \quad u(0)=u_0 \text { on } \Gamma,
\end{equation}\] where \(u\)
satisfies the Kirchhoff vertex conditions \[\begin{equation}
\label{eq:Kcond}
\left\{\phi\in C(\Gamma)\;\Big|\; \forall v\in V:
\sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\}
\end{equation}\]
If \(f=0\), then the solution is
given by \[\begin{equation}
\label{eq:sol_reprentation}
u(s,t) =
\displaystyle\sum_{j\in\mathbb{N}}e^{-\lambda^{\frac{\alpha}{2}}_jt}\left(u_0,
e_j\right)_{L_2(\Gamma)}e_j(s).
\end{equation}\]
# Function to build a tadpole graph and create a mesh
gets_graph_tadpole <- function(h){
edge1 <- rbind(c(0,0),c(1,0))
theta <- seq(from=-pi,to=pi,length.out = 10000)
edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
edges = list(edge1, edge2)
graph <- metric_graph$new(edges = edges)
graph$build_mesh(h = h)
return(graph)
}
Let \(\Gamma_T =
(\mathcal{V},\mathcal{E})\) characterize the tadpole graph with
\(\mathcal{V}= \{v_1,v_2\}\) and \(\mathcal{E}= \{e_1,e_2\}\) as specified in
Figure \(\ref{Interval.Circle.Tadpole}\)c. The left
edge \(e_1\) has length 1 and the
circular edge \(e_2\) has length 2. As
discussed in Subsection \(\ref{subsec:prelim}\), a point on \(e_1\) is parameterized via \(s=\left(e_1, t\right)\) for \(t \in[0,1]\) and a point on \(e_2\) via \(s=\left(e_2, t\right)\) for \(t\in[0,2]\). One can verify that \(-\Delta_\Gamma\) has eigenvalues \(0,\left\{(i \pi / 2)^2\right\}_{i \in
\mathbb{N}}\) and \(\left\{(i \pi /
2)^2\right\}_{2 i \in \mathbb{N}}\) with corresponding
eigenfunctions \(\phi_0\), \(\left\{\phi_i\right\}_{i \in \mathbb{N}}\),
and \(\left\{\psi_i\right\}_{2 i \in
\mathbb{N}}\) given by \(\phi_0(s)=1 /
\sqrt{3}\) and \[\begin{equation*}
\phi_i(s)=C_{\phi, i}\begin{cases}
-2 \sin (\frac{i\pi}{2}) \cos (\frac{i \pi t}{2}), & s \in
e_1, \\
\sin (i \pi t / 2), & s \in e_2,
\end{cases},
\quad
\psi_i(s)=\frac{\sqrt{3}}{\sqrt{2}} \begin{cases}
(-1)^{i / 2} \cos (\frac{i \pi t}{2}), & s \in e_1, \\
\cos (\frac{i \pi t}{2}), & s \in e_2,
\end{cases},
\end{equation*}\] where \(C_{\phi,
i}=1\) if \(i\) is even and
\(C_{\phi, i}=1 / \sqrt{3}\) otherwise.
Moreover, these functions form an orthonormal basis for \(L_2(\Gamma_T)\).
# Function to compute the eigenfunctions
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2])
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2])
if(k==0){
f.e1 <- rep(1,length(x1))
f.e2 <- rep(1,length(x2))
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2)
f.e2 <- sin(pi*k*x2/2)
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
if((k %% 2)==1){
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
f.e2 <- cos(pi*k*x2/2)
f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1])
f <- list(phi=f1,psi=f2/sqrt(3/2))
}
}
return(f)
}
Implementation of \(u\)
time_step <- 0.01
T_final <- 2
time_seq <- seq(0, T_final, by = time_step)
h <- 0.01
graph <- gets_graph_tadpole(h = h)
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.15 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
# Compute the FEM matrices
graph$compute_fem()
x <- graph$mesh$V[, 1]
y <- graph$mesh$V[, 2]
edge_number <- graph$mesh$VtE[, 1]
pos <- sum(edge_number == 1)+1
order_to_plot <- function(v)return(c(v[1], v[3:pos], v[2], v[(pos+1):length(v)], v[2]))
weights <- graph$mesh$weights
kappa <- 1
alpha <- 0.51
n_finite <- 10
EIGENVAL <- c() # initialize empty vector for eigenvalues
EIGENFUN <- NULL # initialize NULL for eigenfunctions matrix
INDEX <- c()
for (j in 0:n_finite) {
lambda_j <- kappa^2 + (j*pi/2)^2
e_j <- tadpole.eig(j,graph)$phi#/sqrt(1000)
EIGENVAL <- c(EIGENVAL, lambda_j) # append scalar to vector
EIGENFUN <- cbind(EIGENFUN, e_j) # append column to matrix
INDEX <- c(INDEX, j)
if (j>0 && (j %% 2 == 0)) {
lambda_j <- kappa^2 + (j*pi/2)^2
e_j <- tadpole.eig(j,graph)$psi#/sqrt(1000)
EIGENVAL <- c(EIGENVAL, lambda_j) # append scalar to vector
EIGENFUN <- cbind(EIGENFUN, e_j) # append column to matrix
INDEX <- c(INDEX, j+0.1)
}
}
coeff <- sample(1:5, length(INDEX), replace = TRUE)
U_0 <- EIGENFUN %*% coeff
#U_0 <- 10*exp(-((x-1)^2 + (y)^2))
n_finite <- c(10,1000)
n_finite1 <- n_finite[1]
U_true1 <- matrix(NA, nrow = length(x), ncol = length(time_seq))
U_true1[, 1] <- U_0
for (k in 1:(length(time_seq)-1)) {
aux_k <- rep(0, length(x))
for (j in 0:n_finite1) {
decay_j <- exp(-time_seq[k+1]*(kappa^2 + (j*pi/2)^2)^(alpha/2))
e_j <- tadpole.eig(j,graph)$phi
aux_k <- aux_k + decay_j*sum(U_0*e_j*weights)*e_j
if (j>0 && (j %% 2 == 0)) {
e_j <- tadpole.eig(j,graph)$psi
aux_k <- aux_k + decay_j*sum(U_0*e_j*weights)*e_j
}
}
U_true1[, k+1] <- aux_k
}
n_finite2 <- n_finite[2]
U_true2 <- matrix(NA, nrow = length(x), ncol = length(time_seq))
U_true2[, 1] <- U_0
for (k in 1:(length(time_seq)-1)) {
aux_k <- rep(0, length(x))
for (j in 0:n_finite2) {
decay_j <- exp(-time_seq[k+1]*(kappa^2 + (j*pi/2)^2)^(alpha/2))
e_j <- tadpole.eig(j,graph)$phi
aux_k <- aux_k + decay_j*sum(U_0*e_j*weights)*e_j
if (j>0 && (j %% 2 == 0)) {
e_j <- tadpole.eig(j,graph)$psi
aux_k <- aux_k + decay_j*sum(U_0*e_j*weights)*e_j
}
}
U_true2[, k+1] <- aux_k
}
x <- order_to_plot(x)
y <- order_to_plot(y)
max_error_at_each_time <- apply(abs(U_true1 - U_true2)/abs(U_true1), 2, max)
U_true1 <- apply(U_true1, 2, order_to_plot)
U_true2 <- apply(U_true2, 2, order_to_plot)
# Create interactive plot
fig <- plot_ly(x = ~time_seq, y = ~max_error_at_each_time, type = 'scatter', mode = 'lines+markers',
line = list(color = 'red', width = 2),
marker = list(size = 4),
name = "Max Error")
fig <- fig %>% layout(title = "Max Error at Each Time Step",
xaxis = list(title = "Time"),
yaxis = list(title = "Max Error"))
plot_data <- data.frame(
x = rep(x, times = ncol(U_true1)),
y = rep(y, times = ncol(U_true1)),
z_true1 = as.vector(U_true1),
z_true2 = as.vector(U_true2),
frame = rep(time_seq, each = length(x))
)
# Compute axis limits
x_range <- range(x)
y_range <- range(y)
z_range <- range(c(U_true1, U_true2))
# Initial plot setup (first frame only)
p <- plot_ly(plot_data, frame = ~frame) %>%
add_trace(
x = ~x, y = ~y, z = ~z_true1,
type = "scatter3d", mode = "lines",
name = paste0("n_finite = ", n_finite1),
line = list(color = "blue", width = 2)
) %>%
add_trace(
x = ~x, y = ~y, z = ~z_true2,
type = "scatter3d", mode = "lines",
name = paste0("n_finite = ", n_finite2),
line = list(color = "red", width = 2)
) %>%
layout(
scene = list(
xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "Value", range = z_range)
),
updatemenus = list(
list(
type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 100, redraw = TRUE), fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
)
)
),
title = "Time: 0"
)
# Convert to plotly object with frame info
pb <- plotly_build(p)
# Inject custom titles into each frame
for (i in seq_along(pb$x$frames)) {
t <- time_seq[i]
err <- signif(max_error_at_each_time[i], 4)
pb$x$frames[[i]]$layout <- list(title = paste0("Time: ", t, " | Max Error: ", err))
}
fig # Display the plot
LS0tCnRpdGxlOiAiU29sdmluZyBhIHBhcmFib2xpYyBlcXVhdGlvbiIKZGF0ZTogIkNyZWF0ZWQ6IDIwLTA0LTIwMjUuIExhc3QgbW9kaWZpZWQ6IGByIGZvcm1hdChTeXMudGltZSgpLCAnJWQtJW0tJVkuJylgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIG1hdGhqYXg6ICJodHRwczovL2Nkbi5qc2RlbGl2ci5uZXQvbnBtL21hdGhqYXhAMy9lczUvdGV4LW1tbC1jaHRtbC5qcyIKICAgIGhpZ2hsaWdodDogcHlnbWVudHMKICAgIHRoZW1lOiBmbGF0bHkKICAgIGNvZGVfZm9sZGluZzogc2hvdyAjIGNsYXNzLnNvdXJjZSA9ICJmb2xkLWhpZGUiIHRvIGhpZGUgY29kZSBhbmQgYWRkIGEgYnV0dG9uIHRvIHNob3cgaXQKICAgIGRmX3ByaW50OiBwYWdlZAogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6CiAgICAgIGNvbGxhcHNlZDogdHJ1ZQogICAgICBzbW9vdGhfc2Nyb2xsOiB0cnVlCiAgICBudW1iZXJfc2VjdGlvbnM6IGZhbHNlCiAgICBmaWdfY2FwdGlvbjogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQphbHdheXNfYWxsb3dfaHRtbDogdHJ1ZQpiaWJsaW9ncmFwaHk6IAogIC0gcmVmZXJlbmNlcy5iaWIKICAtIGdyYXRlZnVsLXJlZnMuYmliCmhlYWRlci1pbmNsdWRlczoKICAtIFxuZXdjb21tYW5ke1xhcn17XG1hdGhiYntSfX0KICAtIFxuZXdjb21tYW5ke1xsbGF2fVsxXXtcbGVmdFx7IzFccmlnaHRcfX0KICAtIFxuZXdjb21tYW5ke1xwYXJlfVsxXXtcbGVmdCgjMVxyaWdodCl9CiAgLSBcbmV3Y29tbWFuZHtcTmNhbH17XG1hdGhjYWx7Tn19CiAgLSBcbmV3Y29tbWFuZHtcVmNhbH17XG1hdGhjYWx7Vn19CiAgLSBcbmV3Y29tbWFuZHtcRWNhbH17XG1hdGhjYWx7RX19CiAgLSBcbmV3Y29tbWFuZHtcV2NhbH17XG1hdGhjYWx7V319Ci0tLQoKYGBge3IgeGFyaW5nYW5FeHRyYS1jbGlwYm9hcmQsIGVjaG8gPSBGQUxTRX0KaHRtbHRvb2xzOjp0YWdMaXN0KAogIHhhcmluZ2FuRXh0cmE6OnVzZV9jbGlwYm9hcmQoCiAgICBidXR0b25fdGV4dCA9ICI8aSBjbGFzcz1cImZhLXNvbGlkIGZhLWNsaXBib2FyZFwiIHN0eWxlPVwiY29sb3I6ICMwMDAwOEJcIj48L2k+IiwKICAgIHN1Y2Nlc3NfdGV4dCA9ICI8aSBjbGFzcz1cImZhIGZhLWNoZWNrXCIgc3R5bGU9XCJjb2xvcjogIzkwQkU2RFwiPjwvaT4iLAogICAgZXJyb3JfdGV4dCA9ICI8aSBjbGFzcz1cImZhIGZhLXRpbWVzLWNpcmNsZVwiIHN0eWxlPVwiY29sb3I6ICNGOTQxNDRcIj48L2k+IgogICksCiAgcm1hcmtkb3duOjpodG1sX2RlcGVuZGVuY3lfZm9udF9hd2Vzb21lKCkKKQpgYGAKCgpgYGB7Y3NzLCBlY2hvID0gRkFMU0V9CmJvZHkgLm1haW4tY29udGFpbmVyIHsKICBtYXgtd2lkdGg6IDEwMCUgIWltcG9ydGFudDsKICB3aWR0aDogMTAwJSAhaW1wb3J0YW50Owp9CmJvZHkgewogIG1heC13aWR0aDogMTAwJSAhaW1wb3J0YW50Owp9Cgpib2R5LCB0ZCB7CiAgIGZvbnQtc2l6ZTogMTZweDsKfQpjb2RlLnJ7CiAgZm9udC1zaXplOiAxNHB4Owp9CnByZSB7CiAgZm9udC1zaXplOiAxNHB4Cn0KLmN1c3RvbS1ib3ggewogIGJhY2tncm91bmQtY29sb3I6ICNmNWY3ZmE7IC8qIExpZ2h0IGdyZXktYmx1ZSBiYWNrZ3JvdW5kICovCiAgYm9yZGVyLWNvbG9yOiAjZTFlOGVkOyAvKiBMaWdodCBib3JkZXIgY29sb3IgKi8KICBjb2xvcjogIzJjM2U1MDsgLyogRGFyayB0ZXh0IGNvbG9yICovCiAgcGFkZGluZzogMTVweDsgLyogUGFkZGluZyBpbnNpZGUgdGhlIGJveCAqLwogIGJvcmRlci1yYWRpdXM6IDVweDsgLyogUm91bmRlZCBjb3JuZXJzICovCiAgbWFyZ2luLWJvdHRvbTogMjBweDsgLyogU3BhY2luZyBiZWxvdyB0aGUgYm94ICovCn0KLmNhcHRpb24gewogIG1hcmdpbjogYXV0bzsKICB0ZXh0LWFsaWduOiBjZW50ZXI7CiAgbWFyZ2luLWJvdHRvbTogMjBweDsgLyogU3BhY2luZyBiZWxvdyB0aGUgYm94ICovCn0KYGBgCgoKTGV0IHVzIHNldCBzb21lIGdsb2JhbCBvcHRpb25zIGZvciBhbGwgY29kZSBjaHVua3MgaW4gdGhpcyBkb2N1bWVudC4KCgpgYGB7cn0KIyBTZXQgc2VlZCBmb3IgcmVwcm9kdWNpYmlsaXR5CnNldC5zZWVkKDE5ODIpIAojIFNldCBnbG9iYWwgb3B0aW9ucyBmb3IgYWxsIGNvZGUgY2h1bmtzCmtuaXRyOjpvcHRzX2NodW5rJHNldCgKICAjIERpc2FibGUgbWVzc2FnZXMgcHJpbnRlZCBieSBSIGNvZGUgY2h1bmtzCiAgbWVzc2FnZSA9IFRSVUUsICAgIAogICMgRGlzYWJsZSB3YXJuaW5ncyBwcmludGVkIGJ5IFIgY29kZSBjaHVua3MKICB3YXJuaW5nID0gVFJVRSwgICAgCiAgIyBTaG93IFIgY29kZSB3aXRoaW4gY29kZSBjaHVua3MgaW4gb3V0cHV0CiAgZWNobyA9IFRSVUUsICAgICAgICAKICAjIEluY2x1ZGUgYm90aCBSIGNvZGUgYW5kIGl0cyByZXN1bHRzIGluIG91dHB1dAogIGluY2x1ZGUgPSBUUlVFLCAgICAgCiAgIyBFdmFsdWF0ZSBSIGNvZGUgY2h1bmtzCiAgZXZhbCA9IFRSVUUsICAgICAgIAogICMgRW5hYmxlIGNhY2hpbmcgb2YgUiBjb2RlIGNodW5rcyBmb3IgZmFzdGVyIHJlbmRlcmluZwogIGNhY2hlID0gRkFMU0UsICAgICAgCiAgIyBBbGlnbiBmaWd1cmVzIGluIHRoZSBjZW50ZXIgb2YgdGhlIG91dHB1dAogIGZpZy5hbGlnbiA9ICJjZW50ZXIiLAogICMgRW5hYmxlIHJldGluYSBkaXNwbGF5IGZvciBoaWdoLXJlc29sdXRpb24gZmlndXJlcwogIHJldGluYSA9IDIsCiAgIyBTaG93IGVycm9ycyBpbiB0aGUgb3V0cHV0IGluc3RlYWQgb2Ygc3RvcHBpbmcgcmVuZGVyaW5nCiAgZXJyb3IgPSBUUlVFLAogICMgRG8gbm90IGNvbGxhcHNlIGNvZGUgYW5kIG91dHB1dCBpbnRvIGEgc2luZ2xlIGJsb2NrCiAgY29sbGFwc2UgPSBGQUxTRQopCmBgYAoKCgoKYGBge3J9CiMgaW5sYS51cGdyYWRlKHRlc3RpbmcgPSBUUlVFKQojIHJlbW90ZXM6Omluc3RhbGxfZ2l0aHViKCJpbmxhYnJ1LW9yZy9pbmxhYnJ1IiwgcmVmID0gImRldmVsIikKIyByZW1vdGVzOjppbnN0YWxsX2dpdGh1YigiZGF2aWRib2xpbi9yc3BkZSIsIHJlZiA9ICJkZXZlbCIpCiMgcmVtb3Rlczo6aW5zdGFsbF9naXRodWIoImRhdmlkYm9saW4vbWV0cmljZ3JhcGgiLCByZWYgPSAiZGV2ZWwiKQpsaWJyYXJ5KElOTEEpCmxpYnJhcnkoaW5sYWJydSkKbGlicmFyeShyU1BERSkKbGlicmFyeShNZXRyaWNHcmFwaCkKbGlicmFyeShncmF0ZWZ1bCkKCmxpYnJhcnkocGxvdGx5KQpgYGAKCgpXZSB3YW50IHRvIHNvbHZlIHRoZSBmcmFjdGlvbmFsIGRpZmZ1c2lvbiBlcXVhdGlvbgpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTptYWluZXF9CiAgICBccGFydGlhbF90IHUrKFxrYXBwYV4yLVxEZWx0YV9cR2FtbWEpXntcZnJhY3tcYWxwaGF9ezJ9fSB1PWYgXHRleHQgeyBvbiB9IFxHYW1tYSBcdGltZXMoMCwgVCksIFxxdWFkIHUoMCk9dV8wIFx0ZXh0IHsgb24gfSBcR2FtbWEsClxlbmR7ZXF1YXRpb259CndoZXJlICR1JCBzYXRpc2ZpZXMgdGhlIEtpcmNoaG9mZiB2ZXJ0ZXggY29uZGl0aW9ucwpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTpLY29uZH0KICAgIFxsZWZ0XHtccGhpXGluIEMoXEdhbW1hKVw7XEJpZ3xcOyBcZm9yYWxsIHZcaW4gVjogXHN1bV97ZVxpblxtYXRoY2Fse0V9X3Z9XHBhcnRpYWxfZSBccGhpKHYpPTAgXHJpZ2h0XH0KXGVuZHtlcXVhdGlvbn0KCklmICRmPTAkLCB0aGVuIHRoZSBzb2x1dGlvbiBpcyBnaXZlbiBieQpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTpzb2xfcmVwcmVudGF0aW9ufQogICAgICAgIHUocyx0KSA9IFxkaXNwbGF5c3R5bGVcc3VtX3tqXGluXG1hdGhiYntOfX1lXnstXGxhbWJkYV57XGZyYWN7XGFscGhhfXsyfX1fanR9XGxlZnQodV8wLCBlX2pccmlnaHQpX3tMXzIoXEdhbW1hKX1lX2oocykuClxlbmR7ZXF1YXRpb259CgpgYGB7cn0KIyBGdW5jdGlvbiB0byBidWlsZCBhIHRhZHBvbGUgZ3JhcGggYW5kIGNyZWF0ZSBhIG1lc2gKZ2V0c19ncmFwaF90YWRwb2xlIDwtIGZ1bmN0aW9uKGgpewogIGVkZ2UxIDwtIHJiaW5kKGMoMCwwKSxjKDEsMCkpCiAgdGhldGEgPC0gc2VxKGZyb209LXBpLHRvPXBpLGxlbmd0aC5vdXQgPSAxMDAwMCkKICBlZGdlMiA8LSBjYmluZCgxKzEvcGkrY29zKHRoZXRhKS9waSxzaW4odGhldGEpL3BpKQogIGVkZ2VzID0gbGlzdChlZGdlMSwgZWRnZTIpCiAgZ3JhcGggPC0gbWV0cmljX2dyYXBoJG5ldyhlZGdlcyA9IGVkZ2VzKQogIGdyYXBoJGJ1aWxkX21lc2goaCA9IGgpCiAgcmV0dXJuKGdyYXBoKQp9CmBgYAoKTGV0ICRcR2FtbWFfVCA9IChcVmNhbCxcRWNhbCkkIGNoYXJhY3Rlcml6ZSB0aGUgdGFkcG9sZSBncmFwaCB3aXRoICRcVmNhbCA9IFx7dl8xLHZfMlx9JCBhbmQgJFxFY2FsID0gXHtlXzEsZV8yXH0kIGFzIHNwZWNpZmllZCBpbiBGaWd1cmUgXHJlZntJbnRlcnZhbC5DaXJjbGUuVGFkcG9sZX1jLiBUaGUgbGVmdCBlZGdlICRlXzEkIGhhcyBsZW5ndGggMSBhbmQgdGhlIGNpcmN1bGFyIGVkZ2UgJGVfMiQgaGFzIGxlbmd0aCAyLiBBcyBkaXNjdXNzZWQgaW4gU3Vic2VjdGlvbiBccmVme3N1YnNlYzpwcmVsaW19LCBhIHBvaW50IG9uICRlXzEkIGlzIHBhcmFtZXRlcml6ZWQgdmlhICRzPVxsZWZ0KGVfMSwgdFxyaWdodCkkIGZvciAkdCBcaW5bMCwxXSQgYW5kIGEgcG9pbnQgb24gJGVfMiQgdmlhICRzPVxsZWZ0KGVfMiwgdFxyaWdodCkkIGZvciAkdFxpblswLDJdJC4gT25lIGNhbiB2ZXJpZnkgdGhhdCAkLVxEZWx0YV9cR2FtbWEkIGhhcyBlaWdlbnZhbHVlcyAkMCxcbGVmdFx7KGkgXHBpIC8gMileMlxyaWdodFx9X3tpIFxpbiBcbWF0aGJie059fSQgYW5kICRcbGVmdFx7KGkgXHBpIC8gMileMlxyaWdodFx9X3syIGkgXGluIFxtYXRoYmJ7Tn19JCB3aXRoIGNvcnJlc3BvbmRpbmcgZWlnZW5mdW5jdGlvbnMgJFxwaGlfMCQsICRcbGVmdFx7XHBoaV9pXHJpZ2h0XH1fe2kgXGluIFxtYXRoYmJ7Tn19JCwgYW5kICRcbGVmdFx7XHBzaV9pXHJpZ2h0XH1fezIgaSBcaW4gXG1hdGhiYntOfX0kIGdpdmVuIGJ5ICRccGhpXzAocyk9MSAvIFxzcXJ0ezN9JCBhbmQgClxiZWdpbntlcXVhdGlvbip9CiAgICBccGhpX2kocyk9Q197XHBoaSwgaX1cYmVnaW57Y2FzZXN9CiAgICAgICAgLTIgXHNpbiAoXGZyYWN7aVxwaX17Mn0pIFxjb3MgKFxmcmFje2kgXHBpIHR9ezJ9KSwgJiBzIFxpbiBlXzEsIFxcClxzaW4gKGkgXHBpIHQgLyAyKSwgJiBzIFxpbiBlXzIsCiAgICBcZW5ke2Nhc2VzfSwKXHF1YWQgCiAgICBccHNpX2kocyk9XGZyYWN7XHNxcnR7M319e1xzcXJ0ezJ9fSBcYmVnaW57Y2FzZXN9CiAgICAoLTEpXntpIC8gMn0gXGNvcyAoXGZyYWN7aSBccGkgdH17Mn0pLCAmIHMgXGluIGVfMSwgXFwKXGNvcyAoXGZyYWN7aSBccGkgdH17Mn0pLCAmIHMgXGluIGVfMiwKXGVuZHtjYXNlc30sClxlbmR7ZXF1YXRpb24qfQp3aGVyZSAkQ197XHBoaSwgaX09MSQgaWYgJGkkIGlzIGV2ZW4gYW5kICRDX3tccGhpLCBpfT0xIC8gXHNxcnR7M30kIG90aGVyd2lzZS4gTW9yZW92ZXIsIHRoZXNlIGZ1bmN0aW9ucyBmb3JtIGFuIG9ydGhvbm9ybWFsIGJhc2lzIGZvciAkTF8yKFxHYW1tYV9UKSQuCgpgYGB7cn0KIyBGdW5jdGlvbiB0byBjb21wdXRlIHRoZSBlaWdlbmZ1bmN0aW9ucyAKdGFkcG9sZS5laWcgPC0gZnVuY3Rpb24oayxncmFwaCl7CngxIDwtIGMoMCxncmFwaCRnZXRfZWRnZV9sZW5ndGhzKClbMV0qZ3JhcGgkbWVzaCRQdEVbZ3JhcGgkbWVzaCRQdEVbLDFdPT0xLDJdKSAKeDIgPC0gYygwLGdyYXBoJGdldF9lZGdlX2xlbmd0aHMoKVsyXSpncmFwaCRtZXNoJFB0RVtncmFwaCRtZXNoJFB0RVssMV09PTIsMl0pIAoKaWYoaz09MCl7IAogIGYuZTEgPC0gcmVwKDEsbGVuZ3RoKHgxKSkgCiAgZi5lMiA8LSByZXAoMSxsZW5ndGgoeDIpKSAKICBmMSA9IGMoZi5lMVsxXSxmLmUyWzFdLGYuZTFbLTFdLCBmLmUyWy0xXSkgCiAgZiA9IGxpc3QocGhpPWYxL3NxcnQoMykpIAogIAp9IGVsc2UgewogIGYuZTEgPC0gLTIqc2luKHBpKmsqMS8yKSpjb3MocGkqayp4MS8yKSAKICBmLmUyIDwtIHNpbihwaSprKngyLzIpICAgICAgICAgICAgICAgICAgCiAgCiAgZjEgPSBjKGYuZTFbMV0sZi5lMlsxXSxmLmUxWy0xXSwgZi5lMlstMV0pIAogIAogIGlmKChrICUlIDIpPT0xKXsgCiAgICBmID0gbGlzdChwaGk9ZjEvc3FydCgzKSkgCiAgfSBlbHNlIHsgCiAgICBmLmUxIDwtICgtMSlee2svMn0qY29zKHBpKmsqeDEvMikKICAgIGYuZTIgPC0gY29zKHBpKmsqeDIvMikKICAgIGYyID0gYyhmLmUxWzFdLGYuZTJbMV0sZi5lMVstMV0sZi5lMlstMV0pIAogICAgZiA8LSBsaXN0KHBoaT1mMSxwc2k9ZjIvc3FydCgzLzIpKQogIH0KfQoKcmV0dXJuKGYpCn0KYGBgCgpJbXBsZW1lbnRhdGlvbiBvZiAkdSQKCmBgYHtyfQp0aW1lX3N0ZXAgPC0gMC4wMQpUX2ZpbmFsIDwtIDIKdGltZV9zZXEgPC0gc2VxKDAsIFRfZmluYWwsIGJ5ID0gdGltZV9zdGVwKQpoIDwtIDAuMDEKZ3JhcGggPC0gZ2V0c19ncmFwaF90YWRwb2xlKGggPSBoKQojIENvbXB1dGUgdGhlIEZFTSBtYXRyaWNlcwpncmFwaCRjb21wdXRlX2ZlbSgpCnggPC0gZ3JhcGgkbWVzaCRWWywgMV0KeSA8LSBncmFwaCRtZXNoJFZbLCAyXQplZGdlX251bWJlciA8LSBncmFwaCRtZXNoJFZ0RVssIDFdCnBvcyA8LSBzdW0oZWRnZV9udW1iZXIgPT0gMSkrMQpvcmRlcl90b19wbG90IDwtIGZ1bmN0aW9uKHYpcmV0dXJuKGModlsxXSwgdlszOnBvc10sIHZbMl0sIHZbKHBvcysxKTpsZW5ndGgodildLCB2WzJdKSkKd2VpZ2h0cyA8LSBncmFwaCRtZXNoJHdlaWdodHMKYGBgCgoKYGBge3J9CmthcHBhIDwtIDEKYWxwaGEgPC0gMC41MQoKCm5fZmluaXRlIDwtIDEwCkVJR0VOVkFMIDwtIGMoKSAgICAgICAjIGluaXRpYWxpemUgZW1wdHkgdmVjdG9yIGZvciBlaWdlbnZhbHVlcwpFSUdFTkZVTiA8LSBOVUxMICAgICAgICMgaW5pdGlhbGl6ZSBOVUxMIGZvciBlaWdlbmZ1bmN0aW9ucyBtYXRyaXgKSU5ERVggPC0gYygpCgpmb3IgKGogaW4gMDpuX2Zpbml0ZSkgewogICAgbGFtYmRhX2ogPC0ga2FwcGFeMiArIChqKnBpLzIpXjIKICAgIGVfaiA8LSB0YWRwb2xlLmVpZyhqLGdyYXBoKSRwaGkjL3NxcnQoMTAwMCkKICAgIEVJR0VOVkFMIDwtIGMoRUlHRU5WQUwsIGxhbWJkYV9qKSAgICAgICAgICMgYXBwZW5kIHNjYWxhciB0byB2ZWN0b3IKICAgIEVJR0VORlVOIDwtIGNiaW5kKEVJR0VORlVOLCBlX2opICAgICAgICAgICAgIyBhcHBlbmQgY29sdW1uIHRvIG1hdHJpeAogICAgSU5ERVggPC0gYyhJTkRFWCwgaikKICAgIGlmIChqPjAgJiYgKGogJSUgMiA9PSAwKSkgewogICAgICBsYW1iZGFfaiA8LSBrYXBwYV4yICsgKGoqcGkvMileMgogICAgICBlX2ogPC0gdGFkcG9sZS5laWcoaixncmFwaCkkcHNpIy9zcXJ0KDEwMDApCiAgICAgIEVJR0VOVkFMIDwtIGMoRUlHRU5WQUwsIGxhbWJkYV9qKSAgICAgICAgICMgYXBwZW5kIHNjYWxhciB0byB2ZWN0b3IKICAgICAgRUlHRU5GVU4gPC0gY2JpbmQoRUlHRU5GVU4sIGVfaikgICAgICAgICAgICAjIGFwcGVuZCBjb2x1bW4gdG8gbWF0cml4CiAgICAgIElOREVYIDwtIGMoSU5ERVgsIGorMC4xKQogICAgfQp9Cgpjb2VmZiA8LSBzYW1wbGUoMTo1LCBsZW5ndGgoSU5ERVgpLCByZXBsYWNlID0gVFJVRSkKVV8wIDwtIEVJR0VORlVOICUqJSBjb2VmZgojVV8wIDwtIDEwKmV4cCgtKCh4LTEpXjIgKyAoeSleMikpCm5fZmluaXRlIDwtIGMoMTAsMTAwMCkKbl9maW5pdGUxIDwtIG5fZmluaXRlWzFdClVfdHJ1ZTEgPC0gbWF0cml4KE5BLCBucm93ID0gbGVuZ3RoKHgpLCBuY29sID0gbGVuZ3RoKHRpbWVfc2VxKSkKVV90cnVlMVssIDFdIDwtIFVfMApmb3IgKGsgaW4gMToobGVuZ3RoKHRpbWVfc2VxKS0xKSkgewogIGF1eF9rIDwtIHJlcCgwLCBsZW5ndGgoeCkpCiAgZm9yIChqIGluIDA6bl9maW5pdGUxKSB7CiAgICBkZWNheV9qIDwtIGV4cCgtdGltZV9zZXFbaysxXSooa2FwcGFeMiArIChqKnBpLzIpXjIpXihhbHBoYS8yKSkKICAgIGVfaiA8LSB0YWRwb2xlLmVpZyhqLGdyYXBoKSRwaGkKICAgIGF1eF9rIDwtIGF1eF9rICsgZGVjYXlfaipzdW0oVV8wKmVfaip3ZWlnaHRzKSplX2oKICAgIGlmIChqPjAgJiYgKGogJSUgMiA9PSAwKSkgewogICAgICBlX2ogPC0gdGFkcG9sZS5laWcoaixncmFwaCkkcHNpCiAgICAgIGF1eF9rIDwtIGF1eF9rICsgZGVjYXlfaipzdW0oVV8wKmVfaip3ZWlnaHRzKSplX2oKICAgIH0KICB9CiAgVV90cnVlMVssIGsrMV0gPC0gYXV4X2sKfQoKbl9maW5pdGUyIDwtIG5fZmluaXRlWzJdClVfdHJ1ZTIgPC0gbWF0cml4KE5BLCBucm93ID0gbGVuZ3RoKHgpLCBuY29sID0gbGVuZ3RoKHRpbWVfc2VxKSkKVV90cnVlMlssIDFdIDwtIFVfMApmb3IgKGsgaW4gMToobGVuZ3RoKHRpbWVfc2VxKS0xKSkgewogIGF1eF9rIDwtIHJlcCgwLCBsZW5ndGgoeCkpCiAgZm9yIChqIGluIDA6bl9maW5pdGUyKSB7CiAgICBkZWNheV9qIDwtIGV4cCgtdGltZV9zZXFbaysxXSooa2FwcGFeMiArIChqKnBpLzIpXjIpXihhbHBoYS8yKSkKICAgIGVfaiA8LSB0YWRwb2xlLmVpZyhqLGdyYXBoKSRwaGkKICAgIGF1eF9rIDwtIGF1eF9rICsgZGVjYXlfaipzdW0oVV8wKmVfaip3ZWlnaHRzKSplX2oKICAgIGlmIChqPjAgJiYgKGogJSUgMiA9PSAwKSkgewogICAgICBlX2ogPC0gdGFkcG9sZS5laWcoaixncmFwaCkkcHNpCiAgICAgIGF1eF9rIDwtIGF1eF9rICsgZGVjYXlfaipzdW0oVV8wKmVfaip3ZWlnaHRzKSplX2oKICAgIH0KICB9CiAgVV90cnVlMlssIGsrMV0gPC0gYXV4X2sKfQpgYGAKCgoKYGBge3J9CnggPC0gb3JkZXJfdG9fcGxvdCh4KQp5IDwtIG9yZGVyX3RvX3Bsb3QoeSkKbWF4X2Vycm9yX2F0X2VhY2hfdGltZSA8LSBhcHBseShhYnMoVV90cnVlMSAtIFVfdHJ1ZTIpL2FicyhVX3RydWUxKSwgMiwgbWF4KQoKVV90cnVlMSA8LSBhcHBseShVX3RydWUxLCAyLCBvcmRlcl90b19wbG90KQpVX3RydWUyIDwtIGFwcGx5KFVfdHJ1ZTIsIDIsIG9yZGVyX3RvX3Bsb3QpCgojIENyZWF0ZSBpbnRlcmFjdGl2ZSBwbG90CmZpZyA8LSBwbG90X2x5KHggPSB+dGltZV9zZXEsIHkgPSB+bWF4X2Vycm9yX2F0X2VhY2hfdGltZSwgdHlwZSA9ICdzY2F0dGVyJywgbW9kZSA9ICdsaW5lcyttYXJrZXJzJywKICAgICAgICAgICAgICAgbGluZSA9IGxpc3QoY29sb3IgPSAncmVkJywgd2lkdGggPSAyKSwKICAgICAgICAgICAgICAgbWFya2VyID0gbGlzdChzaXplID0gNCksCiAgICAgICAgICAgICAgIG5hbWUgPSAiTWF4IEVycm9yIikKCmZpZyA8LSBmaWcgJT4lIGxheW91dCh0aXRsZSA9ICJNYXggRXJyb3IgYXQgRWFjaCBUaW1lIFN0ZXAiLAogICAgICAgICAgICAgICAgICAgICAgeGF4aXMgPSBsaXN0KHRpdGxlID0gIlRpbWUiKSwKICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICJNYXggRXJyb3IiKSkKCgoKcGxvdF9kYXRhIDwtIGRhdGEuZnJhbWUoCiAgeCA9IHJlcCh4LCB0aW1lcyA9IG5jb2woVV90cnVlMSkpLAogIHkgPSByZXAoeSwgdGltZXMgPSBuY29sKFVfdHJ1ZTEpKSwKICB6X3RydWUxID0gYXMudmVjdG9yKFVfdHJ1ZTEpLAogIHpfdHJ1ZTIgPSBhcy52ZWN0b3IoVV90cnVlMiksCiAgZnJhbWUgPSByZXAodGltZV9zZXEsIGVhY2ggPSBsZW5ndGgoeCkpCikKCiMgQ29tcHV0ZSBheGlzIGxpbWl0cwp4X3JhbmdlIDwtIHJhbmdlKHgpCnlfcmFuZ2UgPC0gcmFuZ2UoeSkKel9yYW5nZSA8LSByYW5nZShjKFVfdHJ1ZTEsIFVfdHJ1ZTIpKQoKIyBJbml0aWFsIHBsb3Qgc2V0dXAgKGZpcnN0IGZyYW1lIG9ubHkpCnAgPC0gcGxvdF9seShwbG90X2RhdGEsIGZyYW1lID0gfmZyYW1lKSAlPiUKICBhZGRfdHJhY2UoCiAgICB4ID0gfngsIHkgPSB+eSwgeiA9IH56X3RydWUxLAogICAgdHlwZSA9ICJzY2F0dGVyM2QiLCBtb2RlID0gImxpbmVzIiwKICAgIG5hbWUgPSBwYXN0ZTAoIm5fZmluaXRlID0gIiwgbl9maW5pdGUxKSwKICAgIGxpbmUgPSBsaXN0KGNvbG9yID0gImJsdWUiLCB3aWR0aCA9IDIpCiAgKSAlPiUKICBhZGRfdHJhY2UoCiAgICB4ID0gfngsIHkgPSB+eSwgeiA9IH56X3RydWUyLAogICAgdHlwZSA9ICJzY2F0dGVyM2QiLCBtb2RlID0gImxpbmVzIiwKICAgIG5hbWUgPSBwYXN0ZTAoIm5fZmluaXRlID0gIiwgbl9maW5pdGUyKSwKICAgIGxpbmUgPSBsaXN0KGNvbG9yID0gInJlZCIsIHdpZHRoID0gMikKICApICU+JQogIGxheW91dCgKICAgIHNjZW5lID0gbGlzdCgKICAgICAgeGF4aXMgPSBsaXN0KHRpdGxlID0gIngiLCByYW5nZSA9IHhfcmFuZ2UpLAogICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAieSIsIHJhbmdlID0geV9yYW5nZSksCiAgICAgIHpheGlzID0gbGlzdCh0aXRsZSA9ICJWYWx1ZSIsIHJhbmdlID0gel9yYW5nZSkKICAgICksCiAgICB1cGRhdGVtZW51cyA9IGxpc3QoCiAgICAgIGxpc3QoCiAgICAgICAgdHlwZSA9ICJidXR0b25zIiwgc2hvd2FjdGl2ZSA9IEZBTFNFLAogICAgICAgIGJ1dHRvbnMgPSBsaXN0KAogICAgICAgICAgbGlzdChsYWJlbCA9ICJQbGF5IiwgbWV0aG9kID0gImFuaW1hdGUiLAogICAgICAgICAgICAgICBhcmdzID0gbGlzdChOVUxMLCBsaXN0KGZyYW1lID0gbGlzdChkdXJhdGlvbiA9IDEwMCwgcmVkcmF3ID0gVFJVRSksIGZyb21jdXJyZW50ID0gVFJVRSkpKSwKICAgICAgICAgIGxpc3QobGFiZWwgPSAiUGF1c2UiLCBtZXRob2QgPSAiYW5pbWF0ZSIsCiAgICAgICAgICAgICAgIGFyZ3MgPSBsaXN0KE5VTEwsIGxpc3QobW9kZSA9ICJpbW1lZGlhdGUiLCBmcmFtZSA9IGxpc3QoZHVyYXRpb24gPSAwKSwgcmVkcmF3ID0gRkFMU0UpKSkKICAgICAgICApCiAgICAgICkKICAgICksCiAgICB0aXRsZSA9ICJUaW1lOiAwIgogICkKCiMgQ29udmVydCB0byBwbG90bHkgb2JqZWN0IHdpdGggZnJhbWUgaW5mbwpwYiA8LSBwbG90bHlfYnVpbGQocCkKCiMgSW5qZWN0IGN1c3RvbSB0aXRsZXMgaW50byBlYWNoIGZyYW1lCmZvciAoaSBpbiBzZXFfYWxvbmcocGIkeCRmcmFtZXMpKSB7CiAgdCA8LSB0aW1lX3NlcVtpXQogIGVyciA8LSBzaWduaWYobWF4X2Vycm9yX2F0X2VhY2hfdGltZVtpXSwgNCkKICBwYiR4JGZyYW1lc1tbaV1dJGxheW91dCA8LSBsaXN0KHRpdGxlID0gcGFzdGUwKCJUaW1lOiAiLCB0LCAiIHwgTWF4IEVycm9yOiAiLCBlcnIpKQp9CgoKZmlnICAjIERpc3BsYXkgdGhlIHBsb3QKcGIKCmBgYAoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCg==